function x = JacobiSolver(lhs, rhs, numIter)
    x = zeros(size(lhs, 2), 1);
    D = spdiags(spdiags(lhs, 0), 0, size(lhs, 1), size(lhs, 2));
    LU = lhs - D;
       
    for i = 1 : numIter        
        x = D \ (rhs - LU * x);
        residual = rhs - lhs * x;
        mse = residual' * residual;
       % fprintf(1, 'Iteration %d residual %.3f\n', i, mse);        
    end
end